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Aims. We present a new continuum 3D radiative transfer code, MCFOST, based on a Monte-Carlo method. MCFOST 
can be used to calculate (i) monochromatic images in scattered light and/or thermal emission, (ii) polarisation maps, (iii) 
interferometric visibilities, (iv) spectral energy distributions and (v) dust temperature distributions of protoplanetary disks. 
Methods. Several improvements to the standard Monte Carlo method are implemented in MCFOST to increase efficiency and 
f"**) reduce convergence time, including wavelength distribution adjustments, mean intensity calculations and an adaptive sampling 
\^ , of the radiation field. The reliability and efficiency of the code are tested against a previously defined benchmark, using a 2D 
disk configuration. No significant difference (no more than 10%, and generally much less) is found between the temperatures 
and SEDs calculated by MCFOST and by other codes included in the benchmark. 

Results. A study of the lowest disk mass detectable by Spitzer, around young stars, is presented and the colours of "represen- 
tative" parametric disks are compared to recent IRAC and MIPS Spitzer colours of solar-like young stars located in nearby 
star forming regions. 



Oh' 



o 



Key words, radiative transfer - stars; circumstellar matter - methods : numerical - polarisation - scattering 



^ 1. Introduction young stars. The sensitivity and wavelength range cov- 

ered by new instruments is increasing steadily and the 

Signatures for the presence of dust are found nearly ev- ^^^^ far-infrared ranges are now being explored effi- 

erywhere m astrophysics. In the context of star and planet ^^^^^^^ g^-^^^^ g^^^^ Telescope, and new facilities 

formation, dust is abundant m molecular clouds and m the j^j^^ Herschel and ALMA will soon complete the coverage, 
circumstellar environments of a large fraction of stellar ob- 

jects in the early stages of their evolution. In the circum- With this unprecedented wealth of data, from optical 

, 1, J. 1 • r 4. u 1 J- to radio, fine studies of the dust content and evolution 

stellar disks encircling these young stars, where planets r , , , 

,1 , , ; n 4.U • J- 1 u 4. J 4- J of disks become possible and powerful radiative transfer 

are thought to form, the interplay between dust and gas . ^ , , , r ,i , i , 

r i • i (RT) codes are needed to fully exploit the data. In this 

IS of paramount importance. ^ ' , ^ 

, , 11, . rr. . 1 1 1 paper we describe such a code, MCFOST. 
At short wavelengths, dust grains efhciently absorb, 

scatter, and polarise the starlight. How much radiation is That code was used extensively to produce syn- 

scattered and absorbed is a function of both the geometry t^^etic images of the scattered light from disks around 

of the disk and the properties of the dust. In turns, the young stars. Examples include the circumbinary ring of 

amount of absorbed radiation sets the temperature of the Tau | |McCabe et alj | 2002t Puchene et al. || 2004( ), the 

dust (and gas) and defines the amount of radiation that large silhouette disk associated with IRAS 04158+2805 

is re-emitted at longer, thermal, wavelengths. (Menard et al. 200^ and an analysis of the circular po- 

rrii 1 , 1 1 1 ., J • if larisation in GSS 30 iChrvsostomou et al.lll997l ). 

Ihe last decade has witnessed the improvement ot ' ' " ^ 

imaging capabilities with the advent of potent instruments In §2 below, we briefly describe MCFOST. In see- 
in the optical and near-infrared, and large miUimeter inter- tion §3, tests and vaHdation of MCFOST are presented, 
ferometers, providing detailed views of the disks around Two examples of applications are presented in §4. Firstly, 

a study is presented of the minimum mass of disks de- 

Send offprint requests to: C. Pinte tectable by Spitzer around young solar-like stars, young 

e-mail: christophe.pinte@obs.ujf-grenoble.fr low-mass stars, and young brown dwarfs. Secondly, the 
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colours of parametric disks are compared to the Spitzer 
colours, for both IRAC and MIPS, of samples of T Tauri 
stars located in nearby star forming regions. 

2. Description of the numerical code 

MCFOST is a 3D continuum radiative transfer code based 
on the Monte-Car lo method. It was originally developed 
bv iMengirdl l)l989 ^ to model the scattered light in dense 
dusty media (including linear and circular polarisations). 
In this paper, we present an extended version of the orig- 
inal code that includes dust heating and continuum ther- 
mal re-emission. 

2.1. Geometry of the computation grid 

MCFOST uses a spatial grid that is defined in cylindrical 
coordinates, well adapted to the geometry of circumstellar 
disks. Nr logarithmically-spaced radial grid cells and 
linearly-spaced vertical grid cells are used. The size of the 
vertical cells follows the flaring of the disk, with a cutoff at 
10 times the scale height for each radius. MCFOST allows 
for the density distribution to be defined arbitrarily, in 
3 dimensions, with the limitation that within each cell, 
quantities are held constant. The density distribution of 
the media under consideration can be described either by 
using a parametric model (as in § I3.1|l or by including any 
ad hoc density table, calculated by other means. In the 
following, the term "local" refers to the properties of a 
single cell. 



2.4. Ttie radiative transfer sctieme in MCFOST 

The Monte Carlo method allows to follow individual 
monochromatic "photon packets" that propagate through 
the circumstellar environment until they exit the compu- 
tation grid. The propagation process is governed by scat- 
tering, absorption and re-emission events that are con- 
trolled by the optical properties of the medium (opacity, 
albedo, scattering phase function, etc..) and by the tem- 
perature distribution. Upon leaving the circumstellar en- 
vironment, "photon packets" are used to build an SED 
and/or synthetic images, one image per wavelength. 

2.4.1. The emission of "photon packets" 

In MCFOST, two general radiation sources are consid- 
ered : photospheric emission and circumstellar dust ther- 
mal emission. 

Any number of stars can be considered, at any po- 
sition in the calculation grid. The photospheric emission 
can be represented by (i) a sphere radiating uniformly 
or (ii) a limb-darkened sphere or (iii) a point-like source, 
whichever is relevant for the problem under consideration. 
The spectrum of photospheric photons can follow a black- 
body or any observed or calculated spectrum. Hot and/or 
cold spots can be added to the photospheres. 

For the dust thermal emission, the density, tempera- 
ture and opacities are considered constant within each grid 
cell. Thermal photons emission locations arc uniformly 
distributed within each cell and the emission is assumed 
to be isotropic. 



2.2. Position-dependent dust distributions 

An explicit spatial dependence of the size (and/or compo- 
sition) distribution f{a,r) is implemented in MCFOST. 
The dust properties are therefore defined locally, i.e., each 
cell of the disk may contain its own independent dust pop- 
ulation. This allows for example to model dust settling to- 
wards the disk midplane, variations of the chemical com- 
position of the dust, and/or an increase of ice mantles 
from the inner, hot regions to the outer, cold edge of the 
disk. 

2.3. Optical dust properties 

The dust optical properties are currently computed with 
the Mie theory, i.e., grains are spherical and homogeneous. 
The optical properties at a given position in the disk are 
derived in accordance with the local size and composition 
distributions f{a,r). The extinction and scattering opac- 
ities are given by 

«°^*/--(A, r) = T"""" 7Ta'Q,^,/UX, a)/(a, r) da (1) 

where (5sca(A,a) and Qoxt(A,a) are respectively the scat- 
tering and extinction cross sections of a grain of size a at 
a wavelength A. 



2.4.2. Distance between interactions 

It is natural within a Monte Carlo scheme to estimate the 
distance a photon packet "travels" between two interac- 
tions by means of optical depth (directly related to the 
density) rather than by physical, linear distance. From a 
site of interaction, i.e., scattering and emission, the optical 
depth T to the next site is randomly chosen following the 
probability distribution p{t) = exp(— r). 

The distance / is computed by integrating the infinites- 
imal optical depth k^^^ {X, r) p{r) ds until the following 
equality is verified : 

rx= /^^''*(A,r)p(r)ds . (2) 
Jo 

Once the position of interaction r is determined, the 
probability that the interaction is a scattering event rather 
than an absorption event is estimated with the local 
albedo : 

_ , _ laZT ^"'Qsca(A, a)/(a, r) da 

Jcimin ™ 'yoxt(A, a)/(a,r)da 

2.4.3. Scattering 

MCFOST includes a complete treatment of linear and cir- 
cular polarisations by using the Stokes formalism. The 
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state of a light packet is described by its Stokes vector, 
with I representing the intensity, Q and U the linear po- 
larisation, and V the circular polarisation. The interaction 
of a photon with a dust particle is described by a 4 x 4 
matrix, the Mueller matrix, S, also known as the scat- 
tering matrix. Since the matrix is defined with respect to 
the scattering plane, appropriate rotation matrices are re- 
quired to transform the Stokes parameters between multi- 
ple scattering events. In the simplifying case of Mie scat- 
tering , the matrix becomes block-diagonal with only 4 
different non-zero elements, see Eq. 0] 



( '\ 

Q 

u 



( S\\ Si2 

Sl2 Sii 










'S'33 





534 



\ 


( '\ 

Q 




u 


/ 





(4) 



The direction of scattering is defined by two angles, 
the scattering angle 9 and the azimuth angle (j) in spherical 
coordinates. Each individual element Sij of the matrix is 
a function of the scattering angle and of the wavelength, 
Sij = Sij {9, A). The scattering angle 9 is randomly chosen 
from the pretabulated cumulative distribution function : 



^ j;;Sn{9') sinW 
^ ' Jo Sii{9') sm9'd9' 



(5) 



For linearly unpolarised incoming light packets, i.e., 
Q,U = 0, the distribution of azimuthal angle is 
isotropic. For light packets with a non-zero linear polari- 
sation, P — \/ + U'^ 1 the azimuthal angle is defined 
relative to its direction of polarisation and determined by 
means of the following cumulative distribution function : 

PU^ 1 (. S,m-Si2{9) sin(20) \ 



where 9 was previously chosen from Eg . 1^ l|Solcll 19891) . 
The Mueller matrix used in Eqs. 01 El and El can be 

defined in two different ways : one Mueller matrix per 
grain size or one mean Mueller matrix per grid cell. In the 
latter case, the local Mueller matrix 5(A, r) is defined by : 



S(A,r) = / S(A,a)/(a,r)da 



(7) 



where ^(A, a) is the Mueller matrix of a grain of size a at 
a wavelength A. 

In the former case, the grain size must be chosen ex- 
plicitly for each event following the probability law : 



p{a)Aa 



7ra^Qsca(a)/(Q,T-)da 
J^2^^ 7ra2Qsca(a)/(a, r) da 



(8) 



The two methods are strictly equivalent. They can be used 
alternatively to minimize either the memory space or the 
computation time required. 



2.4.4. Absorption and radiative equilibrium 

The temperature of the dust particles is determined by as- 
suming radiative equilibrium in the whole model volume, 
and by assuming that the dust opacities do not depend on 
temperature. 

The dust thermal balance should take into account the 
thermal coupling between gas and dust. In high density 
regions, e.g., close to the disk midplanc, the coupling is 
very good and the gas temperature is expected to be close 
to the dust temperature. At the surface layers of the disk, 
on the other hand, densities become much lower and gas- 
dust thermal exchanges are reduced. 

Two extreme assumptions are implemented in 
MCFOST : either the gas-dust thermal exchange is perfect 
and gas and dust are in local thermodynamic equilibrium 
(LTE) throughout the disk, or there is no thermal cou- 
pling between gas and dust at all. In the latter case, each 
grain size has its own, different temperature, despite being 
in equilibrium with the same radiation field. 

Under the assumption of LTE, all dust grains have the 
same temperature, which is also equal to the gas temper- 
ature. In this case, the (unique) temperature of the cell is 
given by the radiative equilibrium equation : 



4^ / Kf'^{\)Bx[T,)A\^Vf^ 
'0 



(9) 



where r^'^*' is the energy absorption rate. If only passive 
heating is considered, this rate can be written from the 
mean intensity, leading to : 



47T Kf%\)Bx{Ti)d\ = ATi nf' 
Jo JO 



(A)JAdA. (10) 



Any extra internal source of energy (viscous heating for 
instance) can be easily included by adding the correspond- 
ing term on the right hand side of Eq. 1101 

Each time a packet 7 of a given wavelength A travels 
through a cell, the quantity A/^, the distance travelled by 
the packet in the cell, is computed. The mean intensity Jx 
in the cell is then derived following Lucv. l^flQQ.'l : 



AttV^ ^ ' AttV, ^ 

7 



(11) 



where Vi is the volume of cell i and e the individual lu- 
minosity of a packet. The radiative equilibrium equation 
can then be written as : 



\ ,^abs 



(12) 



This method is known as the mean intensity method. 

When gas and dust, as well as dust grains of different 
sizes, are not thermally coupled, there is a different tem- 
perature for each grain size and the radiative equilibrium 
must be written independently for each of them : 



47r 



abs 



(13) 
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where Kf^^{X,a), Ti{a) and Tf^^{a) are the opacity, tem- 
perature and energy absorption rate of the dust grains of 
size a in cell i. 

The calculation of Kf''{X)Bx{T,) dX = aT^Kp/ir 
{kp is the Planck mean opacity) is very time consum- 
ing and we pretabulate these values at Nt = 1 000 loga- 
rithmically spaced temperatures ranging from IK to the 
sublimation temperature, which is of order 1 500 K. The 
temperature Ti of each cell is obtained by interpolation 
between the pretabulated tem peratures. 

To speed up calculations. iBiorkman &: WoodI l)200l|) 
introduced the concepts of immediate reemission and tem- 
perature correction. They are implemented in MCFOST. 
In essence, when a packet is absorbed, it is immediately 
re-emitted and its wavelength chosen by taking into ac- 
count the temperature correction given by the following 
probability distribution : 



px dX cx nf^iX) 



fdBx 



dX 



(14) 



\dT yrj.^ 

iBaes et all ()2005l) suggested that this wavelength dis- 
tribution adjustment method may not be entirely safe 
when used in combination with the mean intensity writing 
of the absorption rate described above. Instead, MCFOST 
determines the cell temperatures Ti used in Eq.^Jby solv- 



ing : 



Kf%X)Bx{T,)dX 



AirViN. 



-N. 



7 abs,i 



(15) 
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where iV-y abs is the current number of absorbed photons 
in cell i. Eq. ^] is used only at the end, when all pho- 
tons have been propagated, to improve the accuracy of 
the temperature determination. Using Eq. 1151 instead of 
Eq. 1121 when computing temperature for Eq. 1141 is in no 
way a limitation because, at this stage of the computation, 
the aim is to compute the radiation field and not the whole 
temperature structure. Of course, the temperature in op- 
tically thin cells will not be precisely determined because 
little absorption occurs in these cells. On the other hand, 
as a consequence, little reemission will also occur and the 
contribution of these reemissi on events will rema. i n very 
small. By construction of the IBiorkman fc WoodI l)200l[l 
method, temperature is best converged in the cells that 
contribute most to the radiation field. So, when all pho- 
tons have been propagated, temperature corrections are 
no longer needed and t he temperat ure of all cells can be 
computed following the ILucvI |)1999(1 method, providing a 
much higher final accuracy. 

2.4.5. SED calculation: sampling of the radiation field 

To improve its computational efficiency when producing 
SEDs, MCFOST uses different samplings of the radiation 
field, i.e., different strategies to set the energy of a pho- 
ton packet. The energy of a packet, hence the number of 
photons it contains, is chosen and optimised depending on 
the goals of the calculation. On the one hand, the conver- 
gence of the temperature distribution is optimized when 



all photon packets have the same energy, independently of 
their wavelengths. This procedure insures that more pho- 
ton packets are emitted in the more luminous bins. On the 
other hand, the computation of SEDs itself is more effi- 
cient when it is the number of photon packets that is held 
constant for all wavelength bins instead. In that case it is 
the energy of the packets that is wavelength dependent. 
Thus, MCFOST is made more efficient by computing the 
temperature and SEDs with a two-step process: 

— Step 1 is the temperature determination. Photon 
packets are generated and calculated one by one at the 
stellar photosphere and propagated until they exit the 
computation grid, i.e., until they exit the circumstellar 
medium. Upon scattering, the propagation vector of a 
packet is modified, but not its wavelength. Upon ab- 
sorption however, packets are immediately re-emitted, 
in situ and isotropically, but at a different wavelength 
calculated according to the temperature of the grid 
cell. For this re-emission process, the concept of imme- 
diate reemission and the a ssociated temperature cor - 
rection method proposed bv lBiorkman fc Wood! l|2nnl^ 
are used. In this step, all photon packets have the same 
luminosity 

e = L,/N^,tepi (16) 

where is the bolometric luminosity of the star and 
N-ystepi the number of packets generated, and are ran- 
domly scattered/absorbed within the disk. The c on- 
cept of mean intensity suggested by ILucvI l)l999l) is 
further used to compute the final temperature at the 
end of the step. This allows to efficiently reduce noise 
in the temperature estimation for optically thin cases. 
Step 1 allows for a fast convergence of the tempera- 
ture but is rather time-consuming when used to derive 
SEDs, especially in the low-energy, long wavelength 
regime. To speed up convergence, a different scheme 
(step) is used to compute SEDs. 

— Step 2 computes the SED from the final temperature 
distribution calculated in step 1. In step 2, the number 
of photon packets N^step2 is held constant at all wave- 
lengths. Step 2 therefore maintains a comparable noise 
level in each wavelength bins and significantly reduces 
convergence time by limiting the CPU time spent in 
high luminosity bins and focusing on low luminosity 
bins. In this step, photon packets are always scattered, 
no absorption is allowed. However, at each interaction, 
the Stokes vector is weighted by the probability of scat- 
tering, Pscaj to take into account the energy that would 
have been removed by absorption. This allows all the 
photons to exit the disk and to contribute to the SED, 
although with a reduced weight. Packets at a given 
wavelength A are randomly emitted from the star and 
the disk, by cells «, with the respective probabilities : 

L.{X) 



P* 



and p. 



L*W + E» «^«(A)i.(A) ' 
w,iX) UjX) 
L*W + w,{X)L,{X) 



(17) 



(18) 
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where the luminosities are defined by : 

L,{X)=A7:^RlBxin) , 

and L,(A) = in mi nf^X) Bx{T,) . 



(19) 
(20) 



In order to avoid generating photon packets so deep 
in the (very optically thick) disk that none of them 
would escape with an appreciable energy left, a "dark 
zone" is defined for each wavelength. This dark zone 
includes all cells from which the optical depth is at 
least T\ = 30, in any direction, to get out of the com- 
putationnal grid^. No photon packet is emitted from 
this zone, the emission probability from cells inside 
the dark zone being set to zero by fixing a weight 
Wi{X) = 0. A weight Wi{X) = 1 is given to cells 
above the dark zone limit. Furthermore, packets that 
would enter the dark zone during their random walk 
are killed, because they have no chance to exit the 
model volume with significant energy. The luminosity 
of the N^step2 packets emitted at a given wavelength 
A is determined by the total energy that the star and 
the disk radiates at this wavelength : 

L,{X) + w,iX)L,{X) 



ex 



'ystep2 



(21) 



2.4.6. Intensity and polarisation maps 

A scheme similar to step 2 used above to calculate SEDs is 
deployed in MCFOST to calculate synthetic maps. Maps 
are calculated, one wavelength at a time. Scattering is 
treated explicitly for all wavelengths but again, absorption 
is replaced by the appropriate weighting of the intensity to 
avoid losses of photon packets. Images are produced, in all 
four Stokes parameters, by classifying photon packets that 
exit the calculation grid into inclination bins, themselves 
divided into image pixels. The same process applies to 
scattered light images in the optical and thermal emission 
maps in the far infrared and radio range. The relative con- 
tribution of the photospheric and dust thermal emission 
is dependent on the wavelength. It is for instance criti- 
cal in the mid-infrared and beyond where it is the (scat- 
tered) thermal emission from the inner disk that starts to 
dominate over the photospheric contribution. The correct 
fraction is calculated explicitly through the temperature 
and density of each disk cell and each contribution can be 
mapped separately for comparison. Interferometric visibil- 
ities and phases can be obtained by Fourier transform of 
images. 

3. Tests and validation of MCFOST 

Numerous tests have been performed to check the effi- 
ciency, stability, and accuracy of MCFOST. In this sec- 
tion, we present a comparison of MCFOST w it h the 
benchmark results published by IPascucci et alJ l|2004f) 
(hereafter P04). 

^ The value tx ~ 30 was found to be a good compromise 
between result accuracy and computational time. 
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Fig. 1. Plots of the radial temperature (upper panel) 
and the difference (lower panel) for the most optically 
thick case, ry = 100, arbitrarily using RADICAL as 
the reference. In both panels, the results of MCFOST 
are represented by the thick solid line. Thin solid 
lines are the results from MC3D, dot-dashed lines from 
MCTRANSF, dotted lines from RADMC and dashed fines 
from STEINRAY. In the upper panel, because all curves 
are very similar, MCFOST has been shifted by 2GGK for 
clarity. 



3.1. Benchmark definition 

The geometry tested involves a central point-like source 
radiating as a T=5 800K blackbody with a bolometric lu- 
minosity L — Lq encircled by a disk of well defined geome- 
try and dust content. The disk extends from 1 AU to 1 000 
AU. It includes spherical dust grains made of astronomical 
silicate. Grains have a single size radius of 0.12 /itm and 
a density of 3 . 6 g.cm ~^. The optical data are taken from 
iDraine fc Led (jlQS^) . The disk is flared with a Gaussian 
vertical profile. Radial power-law distributions are used 
for the surface density and the scale height (see P04 for 
more details). 

To compare with P04, the scattering is also forced to 
be isotropic and polarisation is not calculated, i.e, all the 
scattering information is contained in the scattering cross 
section alone, Qsca, acting on the I Stokes parameter only. 

The number of grid cells is set to Nr = 50 and Nz = 20 
in the radial and vertical directions, respectively. We used 
10^ photon packets to calculate the temperature distribu- 
tion (step 1) and 2 10^ photon packets per wavelength to 
generate the SED (step 2). The total run-time is 20 min- 
utes for the most optically thick case on a bi-processor 
(Intel Xeon) computer running at a clock rate of 2.4 GHz. 
The run-time memory space needed is 10 MBytes. 
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T,, = 100 




Table 1. Star parameters 
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Fig. 2. Plots of the difference between the model SEDs 
for the most optically thick case, ry — 100, arbitrarily us- 
ing RADICAL as the reference code. The linestyle-to-code 
pairing is the same as in Fig.^ Three different inclination 
angles are considered: i — 12.5° (upper panel), i = 42.5° 
(middle panel), and i = 77.5° (lower panel). 
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3.3. Energy conservation at high optical depth 

As discussed in section 12.4.51 SEDs can be constructed 
from packets escaping the system in both steps of the com- 
putation (the SED from step 2 being better converged by 
construction). As photon packets are generated indepen- 
dently through different processes in step 1 and 2, our 
adaptive sampling of the radiation field method allows an 
a posteriori check of the accuracy of calculations. 

By integrating over all inclination angles and wave- 
lengths, the emerging flux must correspond in both cases 
to the star luminosity (in case of passive heating) . This is 
true by construction for the SED calculated in step 1. 

Assuming a standard flared disk geometry and ISM 
like grains, and varying (via the disk mass) the V band 
equatorial optical depth from 0.1 to 10^, we find that 
the luminosities calculated through the 2 methods always 
agree to within better than 1%. 



3.2. Results 

The resuhs calculated with MCFOST for all 
the test cases presented in P04 are available at 
'http : //www-laog. obs .ujf-grenoble . f r/~cpinte/mcf o 

Only the most optically thick case (tv — 100, see P04), 
i.e., the most difficult one, is presented in Figs.nand[21 

The temperature profiles in the disk and shape of 
the emerging SEDs calculated with MCFOST are in 
excellent agreement with the calculation presented in 
P04. Differences in the radial temperatures computed by 
MCFOST do not exceed 4% with respect to the other 5 
codes in the optically thin case. For the most optically 
thick case, differences are generally of the order of 2-5% 
(see Fig.^. The results from MCFOST are always within 
the range of those produced by other codes. 

Regarding SEDs, MCFOST produces results compa- 
rable to all other codes to better than 10% for low op- 
tical thickness. Deviations as large as 15% are however 
observed in the most defavorable case, i.e., high tilt and 
optical thickness (either at very short wavelength or in 
the 10 fim silicate feature), but again the results lie in the 
range of the values produced by other codes (see Fig. [21). 

Various tests were further performed to show the in- 
dependence of the MCFOST results with respect to grid 
sampling, inclination sampling, and position of the verti- 
cal density cut-off in the disk. Results are not noticeably 
affected. A complete descriptio n of all numerical tests is 
presented in lPinte et al 1 (1200,5^ . 



4. Example of applications of MCFOST 

In this section we present calculations of the SEDs and in- 
frared photometric colours of young objects surrounded by 
circumstellar disks. The colours calculated by MCFOST 
are compared with photometric observations obtained 
^^th the Spitzer Space Telescope. Calculated infrared 
axes as a function of disk mass are also compared with 



Spitzer' s detection limits, as obtained for ex ample dur- 
ing legacy surveys hke c2d l|Evans et al.ll2003|) . The goals 
are (i) to explore whether a disk with "representative" 
geometrical parameters can match the observed infrared 
colours of solar-like young stars and (ii) also to estimate 
the range of disk masses detectable by Spitzer in such pro- 
grammes. 

4.1. Disk geometry and central object properties 

We assume a disk geometry that is representative of the 
geometries deri yed from images of sev eral disk sources 
fHK Tau B : ISt apclfeld t et a,1.l ; HV Tau C 

il2l 



ISta,T.e1feldt et alFoOa: IRAS 04158+2805 riClauser et a,lJ 
200 6 : TW Hya : iKrist et all 120001 ; LkHa 263 C : 
|Ch auvin et al ?'2002'). The disk geometry is kept constant 
through the various model runs but the disk mass is var- 
ied. Only disks with dust masses smaller than 10~^Mq 
are considered in order to focus on the transition between 
optically thick and optically thin disks, i.e., from class II 
to class III young stellar objects. 

The disk is defined by a flared geometry with a vertical 
Gaussian profile p{r,z) = po{r) exp{—z^/2h{r)^). Power- 
law distributions describe the surface density, S(r) = 
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Fig. 3. SEDs for different dust disk masses and central objects, for a pole-oii disk. The left panel presents the results 
for the T Tauri star and i?in = 0.1 AU, the central panel for the very low mass star and the right panel for the brown 
dwarf (defined in table The line type corresponds to the disk dust mass (in solar masses). The horizontal bars 
indicate Spitzer's detection limits in the four IRAC bands and the 24 /im and 71 /im MIPS bands for the molecular 
cloud mapping programme of the Cores to Disks legacy survey. Additional MIPS limits (light gray bars : 170 mJy 
at 24 ^m and 1 000 mJy at 71 /im) for deeper pointed observations of the programme are also plotted at 24 /im and 
71 /im. 



So (r/ro)" , and the scale height, h{r) — Hq (r/ro)^ , where 
r is the radial coordinate in the equatorial plane, ho = 10 
AU is the scale height at the reference radius ro — 100 
AU, /3 — 1.125, and a = —1. The outer radius is set 
to i?out — 300 AU. Two values of the disk inner radius 
are considered: 0.1 AU and 1 AU. A grid of models is 
calculated from dust disk masses ranging from 10~^^ to 
1O~^M0. Calculations are presented with masses sampled 
logarithmically, by steps of \/T0- 

A grain size distribution dn{a) cx a^'^'^da of spherical 
particles ranging from amin = 0.03 nm. to Omax = 1 mm, 
and optical constants from iMathis fc WhiffenI (^98^, 
namely their model A, are used. 

The central stars are assumed to emit like blackbod- 
ies. Three different central light sources are considered. 
They are chosen to represent a typical T Tauri star, a 
low mass star, and a brown dwarf (see Tab. ^ for pa- 
rameters). The radii and effective temperatures of the T 
Tauri and low-ma ss stars are obtained from the models of 
ISiess et all l|2000l) , assuming an age of 2 Myrs and a metal- 
licity Z = 0.1. For the brown dwarf, a mass of 0.05 Mq 
is selected. The phot ospheric properties are taken from 
IChabrier et"al] (|200(/) and from Baraffe et all l)2002() dusty 
models. These parameters are summarized in Tab. ^ 

All models are calculated for a distance d=140 pc, typ- 
ical of nearby star forming regions. 



legacy survey team Cores to Disks for their molecular 
cloud mapping programme ijEvans et alJl200^ ^. 

At all 4 IRAC bands, the photospheric emission is al- 
ways above the detection limit, for all types of objects. At 
24 fXTa, the T Tauri and low-mass photospheres are de- 
tected. However, the photosphere of the brown dwarf is 
« 3 times lower than the detection limit. This is of course 
dependent on the exact detection limit chosen and on the 
actual mass of the brown dwarf. 

For the case considered here, i.e., detection limits 
for large mapping programmes and substellar object of 
O.OSMq, a flux detection at 24 /im readily implies the pres- 
ence of circumstellar material. A more careful analysis is 
required to confirm a disk excess in the case of more mas- 
sive objects. At 71 /im, the photosphere is detected for 
none of the objects and the presence of circumstellar dust 
can be directly inferred from any detection. 

Reverting the argument, for an observing programme 
with such detection limits, Spitzer is expected to detect, 
at 71 fim, dust disk masses as low as 3 10~^ Mq around 
T Tauri stars, with R-.n = 0.1 AU (Fig. 01 left panel). For 
low mass stars, disk masses greater than 3 10~^ Mq will be 
detected (Fig.|3| central panel), and for a brown dwarf of 
O.O5M0, disk masses larger than 3 lO"*" Mq are required 
for detection (Fig. O right panel) . All results are given in 
terms of dust mass. A gas-to-dust ratio has to be assumed 
to extrapolate to total disk masses. 



4.2. Detectability of disk infrared excesses with Spitzer 

In Fig. 13 SEDs are presented as a function of disk dust 
mass for the three objects to evaluate the minimum de- 
tectable disk mass as a function of the luminosity of the 
central star. For comparison, a set of detection limits 
reached by Spitzer for large mapping programmes is over- 
plotted. Here we take the values quoted by the Spitzer 



4.3. Spitzer colours 

iHartmann et alJ 1 2005f) have used Spitzer/IRAC to deter- 
mine the location of TTS (Class II and HI) in colour- 
colour diagrams. Here we attempt to reproduce the global 



^ See the c2d website for more information 
|http : //peggysue . as .utexas . edu/SIRTF/, 
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[4.5] - [5.8] [5.8] - [8.0] 

Fig. 4. IRAC colour-colour diagrams. Black filled circles corresponds to models with R\n — 0.1 AU and large grains 
(imax — 1 mm), gray triangles to models with i?in = 1.0 AU and large grains, white squares to models with i?in — 0.1 
AU, with interstellar dust (omax = 1 /im). Large symbols show inclinations smaller than 80° and small symbols 
i > 80°. Class II and Class III data points from Hartmann et al. (.2005) were converted to surface densities using 
a Gaussian convolution (tr = 0.07 mag). The solid and dashed contours define the location of class II and class III 
objects; respectively. The levels are [0.1, 0.2, 0.4, 0.8] of the maximum of density. 



location of these sources in the colour-colour diagrams. 
The precise location depends on the model geometry. 
Similar studies of IRAC colours have been presented else- 
where (D'Alessi o et a.\..200Q^ . Here, we present the colours 
of a few parametric models of passive disks to explore the 
range of parameters that reproduce IRAC colours. We also 
compa re the same model w" ith MIPS colours recently pub- 
lished ijPadeett et alJl2006|) . 

Fig. m and El show the modelled IRAC and MIPS 
colours for the previously defined T Tauri (IMq) mod- 
els. Models for central objects with lower masses (not 
displayed for clarity) presents similar colours and, as a 
consequence, the comparison of our results for a unique 
central object inass w ith t he distribution of rn asses of 
iHartmann et all l)2005|) and lPadeett et all l)2006() sources 
is relevant. Additional models with interstellar dust, where 
Omax = 1 and larger inner radius = 1 AU are 
also presented. Three data sets ( i?in — 0.1 AU and 
flmax = 1 mm ; _Rin = 1 AU and a„iax = 1 mm ; i?in = 0.1 
AU and Omax — 1 M™) a-^e then displayed, each one cal- 
culated for 15 different disk masses and for 21 different 
inclinations, linearly sampled from pole-on to edge-on. 
Contour plots showing the observed colo urs of Taurus 
class II and class III sources (fr om Hartm ann et alJboOSi 
in Fig.l^and lPadgett et al.l2006l in Fig.O are surimposed. 



I 
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K - [24] 



Fig. 5. MIPS colour-colour diagram. The symbols are the 
same as in Fig. 01 The solid and dashed c ontours define 
the lo cation of class II and class HI from Padgett et alJ 
l|2006l) . respectively. The Gaussian convolution was per- 
formed with (T = 0.07 mag along the H - K axis and 
(J = 0.35 mag along the K -[24] axis. 



Broadly speaking, models sample the two regions of 
the colour-colour diagram where sources are found : low 
disk mass models corresponding to class HI objects and 
higher disk mass models to class II. 



The slight offsets seen between the lowest mass models 
and the class HI region is due to the colours of the black- 
body spectrum we used. Using a more realistic spectrum 
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would produce, for instance, a [3.6] - [4.5] colour index of 
—0.04 mag instead of 0.1 mag for the blackbody spectrum, 
in better agreement with the observations. 

For increasing disk masses, all IRAC and MIPS colour 
indices increase toward the class II region. The models 
also populate the intermediate region where no object are 
observed. These results can be used to confirm that the 
transition between class II (disk dust mass > 3 10^^ M©) 
and class III (disk dust mass < 10^^") Mq must be fast 
enough so that no or very few objects are detected in 
samples of nearly 50 sources (Hartmann ct al. 2005). 

A clear difference is seen between models depending on 
the disk inner radius : models with i?in = 1 AU produce 
redder [4.5] - [5.8] and [5.8] - [8.0] colours. In particular, 
with Ria = 1 AU the [5.8] - [8.0] colour index becomes 
larger than the observed value by about 0.5 mag for large 
masses, suggesting that i?in — 1 AU is an upper limit for 
most class II disks in Taurus. This suggestion is supported 
by the H - K / K - [24] diagram, where no excess is seen 
in H - K for models with Ri^ — 1 AU, well below the 
average H - K = 0.3 seen in observations (see gray triangles 
in Fig. 01 and Ej) . Such a behaviour is explained by the 
much higher temperature reached in models with i?in = 
0.1 AU, for which disk emission peaks between 2 and 3 fim, 
whereas it peaks around 10 /im in models with = 1.0 
AU. 

Models close to edge-on (i > 80°) tend to produce red- 
der colour at large masses than their less tilted counter- 
parts. Because the central star becomes more attenuated 
by its disk, the contribution from the disk emission pro- 
gressively dominates, leading to rising mid-infrared SEDs. 
The resulting colours are reminiscent of Class I objects, 
and highlight the necessity to determine the geometry and 
inclination, in order to correctly assess the nature of a 
class I source, i.e., the presence or not of a massive en- 
velop e (jKenvon fc HartmannlllQQst IWhite fc Hillenbrand! 
l2004^ 

No significant influence of the grains size is detectable 
for models where the star is seen directly (with i < 80°). 
At high inclinations, on the contrary, models with inter- 
stellar dust produce redder colour indices than models 
with large grains, for which absorption tends to be gray. 

For completely edge-on models, when no direct light 
from the star reaches the observer, short wavelength 
colours : [3.6] - [4.5] and H - K strongly depend on scat- 
tering properties. Models with interstellar dust present al- 
most no excess at short wavelength, most of the light being 
scattered starlight. At longer wavelength, direct emission 
from the disk is seen, producing red indices. For mod- 
els with larger grains, the albedo remains non negligible 
in the infrared, and the contribution of thermal emission 
from the central (and hidden) parts of the disk scattered 
by the outer parts towards the observer is non negligible 
and adds to the stellar scattered light, resulting in redder 
[3.6] - [4.5] and H - K indices than with interstellar dust. 



5. Summary 

We have presented a new continuum 3D radiative transfer 
code, MCFOST. The efficiency and rehability of MCFOST 
was tested considering the benchmark configuration de- 
fined by P04. MCFOST was shown to calculate temper- 
ature distributions and SEDs that are in excellent agree- 
ment with previous results from other codes. 

Sets of models for young solar-like stars and low-stars 
are presented and compared to the Spiteer's detection lim- 
its for the programme of the Cores to Disks legacy survey. 
Minimal disks masses of « 10~^ Mq for a T Tauri star 
and w 10^^ Mq for a brown dwarf are needed for the disk 
to be detectable by Spitzer in the mapping mode. 

IRAC and MIPS colours of Taurus Class II and III 
objects are compared to passively heated models with a 
"representative" disk geometry. The average location of 
the two classes of objects are well reproduced, as well as 
the extreme colours of some of the objects, that may cor- 
respond to highly tilted disks. i?in = 1 AU is found to 
be maximum inner radius required to account for the ob- 
served colours of Class II objects. Further modelling of 
individual sources, combining optical, infrared and mil- 
limeter photometry with images and/or IRS spectroscopy, 
will be needed to better understand individual disks, their 
dust properties and their evolution. 
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